SUBROUTINE insertLowerPoints(npt,maxpts,bankL,bankR,ymax,deltZ,y,z)
IMPLICIT NONE
INTEGER, PARAMETER :: prec=KIND(1.0D0)    ! defines precision of problem, # of decimal points
INTEGER, PARAMETER :: iprec=SELECTED_INT_KIND(9)    ! defines precision of integers
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! purpose: gets boundary condition, stin, at time tt
!! parameters:
INTEGER(iprec),INTENT(INOUT) :: npt
INTEGER(iprec),INTENT(IN) :: maxpts  ! sediment solution type
REAL(prec),INTENT(IN) :: bankL,bankR  ! left and right bank
REAL(prec),INTENT(IN) :: ymax  ! if distance between two flat points is larger than ymax, add a points
REAL(prec),INTENT(IN) :: deltZ  ! if distance between two flat points is larger than ymax, add a points
REAL(prec),INTENT(INOUT),DIMENSION(maxpts) :: y,z  ! cross section points
!! error conditions:    returns iflag /= 0
!! side effects:        none known
!! bugs:                none known
!! author:              Jianchun Victor Huang 12-17-2009, USBR
!! references:
!! notes:
!---------------------- local declarations ------------------------------------
INTEGER(iprec) :: k,nInsert,kk,nptNew
REAL(prec) :: y1,y2,y3,z1,z2,z3,dy,dz
REAL(prec) :: thalweg
REAL(prec),PARAMETER :: SMALL=0.01D0

thalweg=9999999.0
DO k=1,npt
    IF (y(k)>=bankL .and. y(k)<=bankR) THEN
        if (z(k)< thalweg) thalweg=z(k)
    ENDIF
ENDDO

! check 3 points thalweg
k=2
DO while (k<npt)
    y1=y(k-1); y2=y(k); y3=y(k+1)
    z1=z(k-1); z2=z(k); z3=z(k+1)
    IF (abs(y2-y1)>yMax .and. abs(y3-y2)>yMax .and. abs(z2-thalweg)<SMALL .and. abs(z2-thalweg)<SMALL .and. abs(z3-thalweg)<SMALL) THEN
        ! insert points
        nInsert=MAX1(ABS(y2-y1)/dmax1(yMax,0.1),ABS(y3-y2)/dmax1(yMax,0.1))
        nptNew=npt+nInsert*2
        IF (nptNew>maxpts) THEN
            write(*,*) "ERROR in insertPints: too much points"
            write(*,*) "maxpts=",maxpts
            write(*,*) "npt=",nptNew
        ENDIF
        ! shift nInsert points
        DO kk=nptNew,k+1+2*nInsert,-1
            y(kk)=y(kk-2*nInsert)
            z(kk)=z(kk-2*nInsert)
        ENDDO
        ! lower z2 by deltZ
        z2=z2-deltZ
        ! interpolate nInsert points fro k-1 to k
        dy=(y2-y1)/real(nInsert+1)
        dz=(z2-z1)/real(nInsert+1)
        DO kk=1,nInsert
            y(k-1+kk)=y1+real(kk)*dy
            z(k-1+kk)=z1+real(kk)*dz
        ENDDO
        ! middle point
        y(k+nInsert)=y2
        z(k+nInsert)=z2
        ! interpolate nInsert points fro k to k+1
        dy=(y3-y2)/real(nInsert+1)
        dz=(z3-z2)/real(nInsert+1)
        DO kk=1,nInsert
            y(k+nInsert+kk)=y2+real(kk)*dy
            z(k+nInsert+kk)=z2+real(kk)*dz
        ENDDO
        npt=nptNew
        k=k+2*nInsert+1
    ELSE
        ! skip to next point
        k=k+1
    ENDIF
ENDDO

! check 2 points thalweg
k=1
DO while (k<npt)
    y1=y(k); y2=y(k+1)
    z1=z(k); z2=z(k+1)
    IF (abs(y2-y1)>yMax .and. abs(z1-thalweg)<SMALL .and. abs(z2-thalweg)<SMALL) THEN
        ! insert points
        nInsert=floor(ABS(y2-y1)/dmax1(yMax*2.0,0.1))*2+1
        nptNew=npt+nInsert
        IF (nptNew>maxpts) THEN
            write(*,*) "ERROR in insertPints: too much points"
            write(*,*) "maxpts=",maxpts
            write(*,*) "npt=",nptNew
        ENDIF
        ! shift nInsert points
        DO kk=nptNew,k+1+nInsert,-1
            y(kk)=y(kk-nInsert)
            z(kk)=z(kk-nInsert)
        ENDDO
        ! interpolate nInsert points
        dy=(y2-y1)/real(nInsert+1)
        dz=(z2-z1)/real(nInsert+1)
        DO kk=1,nInsert
            y(k+kk)=y1+real(kk)*dy
            z(k+kk)=z1+real(kk)*dz-(1.0-ABS((y(k+kk)-0.5*(y1+y2)))/((y2-y1)*0.5))*deltZ
            !y(k+kk) current point y
            ! 0.5*(y1+y2) middle point y
            ! distance between k and k+1
        ENDDO
        npt=nptNew
        k=k+nInsert+1
    ELSE
        ! skip to next point
        k=k+1
    ENDIF
ENDDO

END SUBROUTINE insertLowerPoints


